				  
$title  CO2_ESTIMATIONS.GMS

* ###################################

* Per capita income, consumption patterns, and CO2 emissions

* Journal of the Association of Environmental and Resource Economists

* Justin Caron and Thibault Fally

* April 2021

*###################################


* THIS FILE COMPUTES ALL CO2 / ENERGY RELATED STATISTICS


$set SLASH \
$if %system.filesys% == UNIX $set SLASH /

* load gravity data from stata or gams ?
* choices : stata, gams
$ if not set gravitydata $set gravitydata gams


$if not set datadir $set datadir "..%SLASH%data%SLASH"
$setglobal datadir %datadir%

* --------------------------------------------
* DEFINE CALIBRATION TYPE

* -- set of sectors - options : all, fdonly
$if not set sectorset $set sectorset all

* -- type of factor intensity - options : avg, cs - if using average, should change definiton of ENDOW
$if not set beta $set beta cs

* HOMOTHETHIC OR NON-HOMOTHETIC?

* in this file, the "content" parameters are computed for H and NH.

* for the simulation results, though, run either H or NH then combine


* preferences: CLM OR CRIE?

* CRIE
$if not set spec $set spec NH
$if not set demandest $set demandest tc


** -- IF CHOSING FITTED DEMAND FITDEMFLAG
* options: yes / no

* NOTE: for the "content" computation, set to no
$if not set fitdem $set  fitdem yes

* --------------------

* only relevant for the simulation results
* declare supply side specification: noresprod, resprod, resprod_06, resprod_15,  resprod_50
*&$if not set prodspec $set  prodspec noresprod
*$if not set prodspec $set  prodspec resprod_06
$if not set prodspec $set  prodspec resprod_15
*$if not set prodspec $set  prodspec resprod_50


* -- declare dataset
$if not set ds $set ds gtap8gas
*$if not set ds $set ds gtap7_all


* LOAD IO coefficients and other things from CO2preparation.gms
* LOAD SIMULATION RESULTS FROM CF_CO2_NH.GMS	

Sets
r country  /
AUS,NZL,CHN,HKG,JPN,KOR,MNG,TWN,KHM,IDN,LAO,MYS,PHL,SGP,THA,VNM,BGD,IND,NPL,PAK,LKA,CAN,
USA,MEX,ARG,BOL,BRA,CHL,COL,ECU,PRY,PER,URY,VEN,CRI,GTM,HND,NIC,PAN,SLV,AUT,BEL,CYP,CZE,
DNK,EST,FIN,FRA,DEU,GRC,HUN,IRL,ITA,LVA,LTU,LUX,MLT,NLD,POL,PRT,SVK,SVN,ESP,SWE,GBR,CHE,
NOR,ALB,BGR,BLR,HRV,ROU,RUS,UKR,KAZ,KGZ,ARM,AZE,GEO,BHR,IRN,KWT,QAT,SAU,TUR,ARE,EGY,MAR,
TUN,CMR,CIV,GHA,NGA,SEN,ETH,KEN,MDG,MWI,MUS,MOZ,TZA,UGA,ZMB,ZWE,BWA,NAM,ZAF,ISR,OMN
/;

alias (r,s);
alias (r,rloop);

Sets
i industry  /
isr,obs,ros,osg,pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,frs,fsh,coa,oil,gas,omn,cmt,omt,vol,mil,pcr,
sgr,ofd,b_t,tex,wap,lea,lum,ppp,p_c,crp,nmm,i_s,nfm,fmp,mvh,otn,ele,ome,omf,ely,wtr,cns,trd,otp,wtp,atp,cmn,ofi
/;

alias (i,j);

Sets
f factor /Land,UnskLab,SkLab,Capital,NatlRes/;

* for aggregation in the decompositions
SET  ii_	 Sectoral classifications /

MAN	Manufactured and Processed Goods
FOS
ELE
SRV
TRN
AGR

/;


SET	mapi(i,*)	Mapping from GTAP sectors to EPPA-10 sectors /
ATP.TRN,
B_T.AGR,
C_B.AGR,
CMN.MAN,
CMT.AGR,
CNS.MAN,
COA.FOS,
CRP.MAN,
OIL.FOS,
CTL.AGR,
*DWE.SRV,
ELE.MAN,
ELY.ELE,
FMP.MAN,
FRS.AGR,
FSH.AGR,
GAS.FOS,
PDR.AGR,
WHT.AGR,
I_S.MAN,
ISR.SRV,
LEA.MAN,
LUM.MAN,
 MIL.AGR,
MVH.MAN,
NFM.MAN,
NMM.MAN,
OAP.AGR,
OBS.SRV,
OCR.AGR,
OFD.AGR,
OFI.SRV,
OME.MAN,
OMF.MAN,
OMN.MAN,
OMT.AGR,
OSD.AGR,
OSG.SRV,
OTN.MAN,
OTP.TRN,
P_C.FOS,
PCR.AGR,
PFB.AGR,
PPP.MAN,
ROS.SRV,
SGR.AGR,
TEX.MAN,
TRD.SRV,
V_F.AGR,
VOL.AGR,
WAP.MAN,
WTP.TRN,
WTR.SRV
/;


* ------------------------------------
* IF IMPORT GRAVITY AND DEMAND ESTIMATES

* gravity
parameter	coeffs
		ex,im,
		cst
		phiest;

* demand
parameter  lambdaest_;

parameter
*fe(i)           previously estimated demand fixed effects
sigma(i)        previously estimated sigmas
scale           scaling parameter for sigma
technology(i,r) export fixed effects from GE
tcost(i,r,s)    fitted transport costs from GE
theta           calibrated using estimates from the literature?
fittedexp	fitted expenditure
incelast
;


* load Demand estiamtes

*$gdxin estimates\estimates_%ds%_consshare_%demandest%_rall.gdx

* NORMAL CRIE CASE
$if "%spec%"=="NH" $gdxin estimates\estimates_%ds%_logweighted_%demandest%_rall.gdx
$if "%spec%"=="H" $gdxin estimates\estimates_%ds%_logweighted_%demandest%_rall.gdx

* CLM CASE
$if "%spec%"=="NH_CLM" $gdxin estimates\estimates_%ds%_logweighted_%demandest%_CLM.GDX
$if "%spec%"=="NH_CLMshift" $gdxin estimates\estimates_%ds%_logweighted_%demandest%_CLM_flex_shifter.GDX
$if "%spec%"=="NH_CLMquad" $gdxin estimates\estimates_%ds%_logweighted_%demandest%_CLM_flex_quadratic.GDX

$load  coeffs, ex, im, cst,  lambdaest_=lambda, scale= sigma_scale fittedexp incelast


incelast("total",i,r) = 0;

* this is the nh scale
sigma(i) = scale("non-homoth") * coeffs("sigma","coeff",i);

theta(i) = coeffs("backed out theta","coeff",i);
* set theta = 4 for sectors for which we have no estimates
theta(i)$(not theta(i)) = 4;

* -----------------------------------------------------------------------
* LOAD INPUT OUTPUT AND OTHER COEFFICIENTS FROM DATAPREPARATION.GMS

*parameter  beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop;

parameter  beta, intersh, co2int, beta_bilat_total, co2coeff,  trade, fd_true, fd, pcexp, incelast, bmkco2, energy, production, pop, energyint, fd_disagg;
parameter incelast_;


* And everything else from DATAPREPARATION.GMS

* for reporting, load country specific (CS) input-output coefficients

* if doing theta8, have to comment this
$gdxin estimates%SLASH%co2_data_%demandest%_cs.gdx

$if "%demandest%"=="theta8" $gdxin estimates%SLASH%co2_data_theta4_cs.gdx

$load  beta, intersh, co2int,  co2coeff, trade, fd_true=fd, pcexp,  bmkco2, energy, production, pop, energyint, fd_disagg incelast_=incelast

incelast("total",i,r) = incelast_("total",i,r) ;


fd("data",i,r) = fd_true(i,r); 

* note: fittedexp comes from the demand.gms file
fd("non-homoth",i,r) = fittedexp("non-homoth",i,r); 
fd("homoth",i,r) = fittedexp("homoth",i,r); 
fd("homoth partial",i,r) = fittedexp("homoth partial",i,r); 


* -----------------------------------------------------------------------
* COMPUTE MORE COEFFICIENTS AND INTERMEDIATE VALUES

* note: may not all be necessary here

parameter
dem_sh(i,r)  share of sector i in total final demand
trade_sh(i,r,s)  share of country r in country s's absorption (Import shares)
prod_sh(i,r,s)  share of country s in country r's sales
fact_sh(f,i,r)   share of sector i in demand for factor f
endow_sh(f,r) share of factor f in total endowment (income)
gosh_m(j,i,r) gosh matrix coefficient (how much j is used in industry i)
gosh_f(j,r) gosh matrix coefficient (how much j is used by final consumers);

parameters production, absorption;

* --------------------------------------------------------------------------------------------------------------------


$if "%fitdem%"=="yes"  fd("data",i,r) = fd("non-homoth",i,r);

*recompute income elasticity -- CRIE

$if "%fitdem%"=="yes" IncElast("direct",i,r) = incelast("DIRECT fitted dem",i,r) ;


dem_sh(i,r) =  fd("data",i,r) / sum(i.local, fd("data",i,r));


* defined as: from r to s:
trade_sh(i,r,s)$sum(r.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(r.local , trade("data",i,r,s));


* simply loaded
*production(i,r) = sum(s, trade("data",i,r,s));
* INSTEAD OF:
*production(i,r) = vom(i,r);

absorption(j,r) = fd("data",j,r) + sum(i,intersh(j,i,r) * production(i,r));


* -----------------------------------------------------------------------
* LOAD SIMULATIONS RESULTS FROM CF_CO2_NH.GMS


positive variables
                chg_income(r) per capita income
                chg_lambda(r) lagrance multiplier
                chg_fd(i,r) final demand
                chg_absorption(i,r) absorption
                chg_production(i,r) production
                chg_Index(i,r) Price index P
                chg_supplier(i,r) Supplier effect S
                chg_fprice(f,r) factor price
*                chg_wage(f_Lab,r) wages

* only relevant for the sectors in res_sect
		chg_resprice(i,r)  price of sector specific resource
               
* intermediate nests to simplify notation (not necessary)
		res_nonres(i,r)   CES
		factorcost(i,r)	  C-D;
	


parameter supply_elast  calibrated supply elasticity,
		income_shock;

*$gdxin 'results/CFresults_co2_%spec%_%demandest%.gdx'

* REVISION:  CHOSE HERE
*$if "%fitdem%"=="yes" $gdxin 'results/CFresults_co2_%spec%_%demandest%_%prodspec%_fitteddem_10pctshock.gdx'

* 'NORMAL' CASE:
$if "%fitdem%"=="yes" $gdxin 'results%SLASH%CFresults_co2_%spec%_%demandest%_%prodspec%_fitteddem.gdx'

* or, if using observed GDP growth:
*parameter income_Shock_r;
*$if "%fitdem%"=="yes" $gdxin 'results/CFresults_co2_%spec%_tc_%prodspec%_fitteddem_obsgdp.gdx'

* CHANGE INCOME_SHOCK BELOW AS WELL

$if "%fitdem%"=="no" $gdxin 'results%SLASH%CFresults_co2_%spec%_%demandest%_%prodspec%_fitteddem.gdx'
$load supply_elast,  chg_income, chg_lambda, chg_fd, chg_absorption, chg_production,  chg_Index, chg_supplier, chg_fprice, chg_resprice 


$load income_shock

parameter                 chg_production_qty(r) production;


parameter share, sharetotal ;


* -----------------------------------------------------------------------
* ENERGY SET DEFINITION

Set	ene_all(i) all energy goods / ely, coa, oil, gas, p_c /,
	ff(i)   Fossil   fuels       /p_c, oil, gas, coa/     ,
	pe(i)	primary energy        /coa,oil,gas/   
	se(i)	secondary energy        /coa,ely,gas, p_c/  
        e(i)    Energy inputs        /coa,p_c,oil,gas,ely/ ;


* redefining income classifications
set lowincome(r), lowincome2;
set middleincome(r), middleincome2;
set highincome(r);

$ontext
* using the world bank guidelines for analytical classification in 2007 (GNI/capita):
 https://datahelpdesk.worldbank.org/knowledgebase/articles/378833-how-are-the-income-group-thresholds-determined


Low income
Lower middle income
Upper middle income
High income

<= 935
936-3,705
3,706-11,455
> 11,455
$offtext

* lowincome
lowincome(r)$(pcexp(r) < 935) = 1;

* lower middle income
lowincome2(r)$(pcexp(r) > 935 and pcexp(r) < 3705) = 1;

* upper middle income
middleincome2(r)$(pcexp(r) > 3705 and pcexp(r) < 11455) = 1;

* aggregating middle income
middleincome(r)$(pcexp(r) > 935 and pcexp(r) < 11455) = 1;

highincome(r)$(pcexp(r) > 11455) = 1;


* ---------------------------------------------------------------------
* CO2 AND ENERGY ACCOUNTING

* can be combined with co2 later on
parameter	beta_secen   matrix defining presence of secondary energy demand
		secen_bilat_sh
		secen_fd_sh
		content  baseline values
		;

set fdtype / data, non-homoth, homoth, "homoth partial"/;

ALIAS(i,i_) ; 
alias(r,r_);

* computing average MRIO intensity

* COPY THE "CONTENT" PARAMTER AND AVGCO2CONTENT

* do reporting for the elements of the following set:

* note, can add secondary energy to this

* can also add nh4, no2, etc if desired
set gas / co2, ghg, secen /;


* replacing everything into the same parameter
co2int("secen","MRIO output",i,r)    =   energyint("secondary","MRIO output",i,r);
co2int("secen","MRIO output",i,"avg") = energyint("secondary","MRIO output",i,"avg") ;
			   
			   
co2int("secen","MRIO cons total",i,r)    =      energyint("secondary","MRIO cons total",i,r);
co2int("secen","MRIO cons total",i,"avg") =   energyint("secondary","MRIO cons total",i,"avg");

co2int("secen","output",i,r) = energyint("secondary","output",i,r) ;
co2int("secen","output",i,"avg") = energyint("secondary","output",i,"avg") ;

* HERE: PRODUCTION

content(gas,"cons - direct", "total fd", r) = sum(i, fd("data",i,r));
content(gas,"cons - direct", "pop", r) = pop(r);
content(gas,"cons - direct", "logpci", r) = log(pcexp(r));
content("co2","cons - direct - avg int", fdtype, r) = sum(i,co2coeff(i,"fd","all") * fd(fdtype,i,r)) / sum(i, fd(fdtype,i,r));
content("co2","cons - direct - avg int", fdtype, "lowincome") = sum((i,lowincome), co2coeff(i,"fd","all") * fd(fdtype,i,lowincome)) / sum((i,lowincome), fd(fdtype,i,lowincome));
content("co2","cons - direct - avg int", fdtype, "lowincome2") = sum((i,lowincome2), co2coeff(i,"fd","all") * fd(fdtype,i,lowincome2)) / sum((i,lowincome2), fd(fdtype,i,lowincome2));
content("co2","cons - direct - avg int", fdtype, "middleincome2") = sum((i,middleincome2), co2coeff(i,"fd","all") * fd(fdtype,i,middleincome2)) / sum((i,middleincome2), fd(fdtype,i,middleincome2));
content("co2","cons - direct - avg int", fdtype, "highincome") = sum((i,highincome), co2coeff(i,"fd","all") * fd(fdtype,i,highincome)) / sum((i,highincome), fd(fdtype,i,highincome));
content("co2","cons - direct - avg int", fdtype, "middleincome") = sum((i,middleincome), co2coeff(i,"fd","all") * fd(fdtype,i,middleincome)) / sum((i,middleincome), fd(fdtype,i,middleincome));
content("co2","cons - direct", fdtype, r_) =  sum(i_,co2coeff(i_,"fd",r_) * fd(fdtype,i_,r_)) / sum(i_, fd(fdtype,i_,r_));
content("co2","cons - direct", fdtype, "lowincome2") = sum((i,lowincome2), co2coeff(i,"fd",lowincome2) * fd(fdtype,i,lowincome2)) / sum((i,lowincome2), fd(fdtype,i,lowincome2));
content("co2","cons - direct", fdtype, "middleincome2") = sum((i,middleincome2), co2coeff(i,"fd",middleincome2) * fd(fdtype,i,middleincome2)) / sum((i,middleincome2), fd(fdtype,i,middleincome2));
content("co2","cons - direct", fdtype, "lowincome") = sum((i,lowincome), co2coeff(i,"fd",lowincome) * fd(fdtype,i,lowincome)) / sum((i,lowincome), fd(fdtype,i,lowincome));
content("co2","cons - direct", fdtype, "middleincome") = sum((i,middleincome), co2coeff(i,"fd",middleincome) * fd(fdtype,i,middleincome)) / sum((i,middleincome), fd(fdtype,i,middleincome));
content("co2","cons - direct", fdtype, "highincome") = sum((i,highincome), co2coeff(i,"fd",highincome) * fd(fdtype,i,highincome)) / sum((i,highincome), fd(fdtype,i,highincome));

* with electricity:

co2coeff(i,"fd-ele",r_) = co2coeff(i,"fd",r_);
co2coeff("ely","fd-ele",r_) = co2int("co2","MRIO cons total","ely",r_);

co2coeff(i,"fd-ele","all") = co2coeff(i,"fd","all");
co2coeff("ely","fd-ele","all") = co2int("co2","MRIO cons total","ely","avg");

* if using fitted demand shares, need to assign values to some countries which have zero:
* with observed demand, shouldn't matter
$if "%fitdem%"=="yes" co2coeff(i,"fd-ele",r_)$(not co2coeff(i,"fd-ele",r_)) = co2coeff(i,"fd-ele","all");

co2coeff(i,"ghg fd-ele",r_) = co2coeff(i,"fd",r_);
co2coeff("ely","ghg fd-ele",r_) = co2int("ghg","MRIO cons total","ely",r_);

co2coeff(i,"ghg fd-ele","all") = co2coeff(i,"fd","all");
co2coeff("ely","ghg fd-ele","all") = co2int("ghg","MRIO cons total","ely","avg");

content("co2","cons - direct ele", fdtype, r) = sum(i,co2coeff(i,"fd-ele",r) * fd(fdtype,i,r)) / sum(i, fd(fdtype,i,r));
content("co2","cons - direct ele - avg int", fdtype, r) = sum(i,co2coeff(i,"fd-ele","all") * fd(fdtype,i,r)) / sum(i, fd(fdtype,i,r));

* using household FD only  -- does not substantially improve the results
content("co2","cons - direct ele - hh only", "data", r) = sum(i,co2coeff(i,"fd-ele",r) *fd_disagg("c",i,r)) / sum(i, fd_disagg("c",i,r));
content("co2","cons - direct ele - hh only - avg int", "data", r) = sum(i,co2coeff(i,"fd-ele","all") * fd_disagg("c",i,r)) / sum(i, fd_disagg("c",i,r));

content("co2","cons - direct ele - avg int", fdtype, "lowincome2") = sum((i,lowincome2), co2coeff(i,"fd-ele","all") * fd(fdtype,i,lowincome2)) / sum((i,lowincome2), fd(fdtype,i,lowincome2));
content("co2","cons - direct ele - avg int", fdtype, "middleincome2") = sum((i,middleincome2), co2coeff(i,"fd-ele","all") * fd(fdtype,i,middleincome2)) / sum((i,middleincome2), fd(fdtype,i,middleincome2));

content("co2","cons - direct ele - avg int", fdtype, "lowincome") = sum((i,lowincome), co2coeff(i,"fd-ele","all") * fd(fdtype,i,lowincome)) / sum((i,lowincome), fd(fdtype,i,lowincome));
content("co2","cons - direct ele - avg int", fdtype, "middleincome") = sum((i,middleincome), co2coeff(i,"fd-ele","all") * fd(fdtype,i,middleincome)) / sum((i,middleincome), fd(fdtype,i,middleincome));
content("co2","cons - direct ele - avg int", fdtype, "highincome") = sum((i,highincome), co2coeff(i,"fd-ele","all") * fd(fdtype,i,highincome)) / sum((i,highincome), fd(fdtype,i,highincome));

content("co2","cons - direct ele", fdtype, "lowincome2") = sum((i,lowincome2), co2coeff(i,"fd-ele",lowincome2) * fd(fdtype,i,lowincome2)) / sum((i,lowincome2), fd(fdtype,i,lowincome2));
content("co2","cons - direct ele", fdtype, "middleincome2") = sum((i,middleincome2), co2coeff(i,"fd-ele",middleincome2) * fd(fdtype,i,middleincome2)) / sum((i,middleincome2), fd(fdtype,i,middleincome2));

content("co2","cons - direct ele", fdtype, "lowincome") = sum((i,lowincome), co2coeff(i,"fd-ele",lowincome) * fd(fdtype,i,lowincome)) / sum((i,lowincome), fd(fdtype,i,lowincome));
content("co2","cons - direct ele", fdtype, "middleincome") = sum((i,middleincome), co2coeff(i,"fd-ele",middleincome) * fd(fdtype,i,middleincome)) / sum((i,middleincome), fd(fdtype,i,middleincome));
content("co2","cons - direct ele", fdtype, "highincome") = sum((i,highincome), co2coeff(i,"fd-ele",highincome) * fd(fdtype,i,highincome)) / sum((i,highincome), fd(fdtype,i,highincome));

*** ONLY FOR DIRECT FOR NOW: DO THE SAME FOR INDIREC (MORE COMPLICATED)
* using "average tech" i.e. region specific prices

*co2coeff("ely","fd avg tech",r) = co2int("co2","MRIO cons","ely","avg");
*content("co2","cons - direct ele - avg int", fdtype, r) = sum(i,co2coeff(i,"fd avg tech",r) * fd(fdtype,i,r)) / sum(i, fd(fdtype,i,r));


* with GHG
*direct = same as co2
content("ghg","cons - direct - avg int", fdtype, r) = sum(i,co2coeff(i,"fd","all") * fd(fdtype,i,r)) / sum(i, fd(fdtype,i,r));
content("ghg","cons - direct", fdtype, r_) =  sum(i_,co2coeff(i_,"fd",r_) * fd(fdtype,i_,r_)) / sum(i_, fd(fdtype,i_,r_));

content("ghg","cons - direct ele", fdtype, r) = sum(i,co2coeff(i,"ghg fd-ele",r) * fd(fdtype,i,r)) / sum(i, fd(fdtype,i,r));
content("ghg","cons - direct ele - avg int", fdtype, r) = sum(i,co2coeff(i,"ghg fd-ele","all") * fd(fdtype,i,r)) / sum(i, fd(fdtype,i,r));


content(gas,"cons - mrio", fdtype, r_) = 	content(gas,"cons - direct", fdtype, r_) + sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,r_) )  / sum((i), fd(fdtype,i,r_)) ;

* by income level:
content(gas,"cons - mrio", fdtype, "lowincome2") = content(gas,"cons - direct", fdtype, "lowincome2") + sum((i,lowincome2), fd(fdtype,i,lowincome2)* co2int(gas,"MRIO cons indirect",i,lowincome2) )  / sum((i,lowincome2), fd(fdtype,i,lowincome2)) ;
content(gas,"cons - mrio", fdtype, "middleincome2") = content(gas,"cons - direct", fdtype, "middleincome2") + sum((i,middleincome2), fd(fdtype,i,middleincome2)* co2int(gas,"MRIO cons indirect",i,middleincome2) )  / sum((i,middleincome2), fd(fdtype,i,middleincome2)) ;
content(gas,"cons - mrio", fdtype, "lowincome") = content(gas,"cons - direct", fdtype, "lowincome") + sum((i,lowincome), fd(fdtype,i,lowincome)* co2int(gas,"MRIO cons indirect",i,lowincome) )  / sum((i,lowincome), fd(fdtype,i,lowincome)) ;
content(gas,"cons - mrio", fdtype, "middleincome") = content(gas,"cons - direct", fdtype, "middleincome") + sum((i,middleincome), fd(fdtype,i,middleincome)* co2int(gas,"MRIO cons indirect",i,middleincome) )  / sum((i,middleincome), fd(fdtype,i,middleincome)) ;
content(gas,"cons - mrio", fdtype, "highincome") = content(gas,"cons - direct", fdtype, "highincome") + sum((i,highincome), fd(fdtype,i,highincome)* co2int(gas,"MRIO cons indirect",i,highincome) )  / sum((i,highincome), fd(fdtype,i,highincome)) ;



* indirect only

content(gas,"cons - mrio indirect", fdtype, r_) = 
	 sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,r_) )
		  / sum((i), fd(fdtype,i,r_)) ;

* indirect no ele
content(gas,"cons - mrio indirect no ele", fdtype, r_) = 
	 sum(i$(not sameas(i,"ely")), fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,r_) )
		  / sum((i), fd(fdtype,i,r_)) ;



content(gas,"cons - mrio indirect - avg int", fdtype, r_) = 
	 sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") )
		  / sum((i), fd(fdtype,i,r_)) ;


content(gas,"cons - mrio indirect no ele - avg int", fdtype, r_) = 
	 sum(i$(not sameas(i,"ely")), fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") )
		  / sum((i), fd(fdtype,i,r_)) ;

* compute total
content(gas,"cons - mrio indirect - avg int", fdtype, "total") = 
	 sum((r_,i), fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") );




* MRIO total at average int
content(gas,"cons - mrio - avg int", fdtype, r_) = 
	content(gas,"cons - direct - avg int", fdtype, r_) +
	 sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") )
		  / sum((i), fd(fdtype,i,r_)) ;

* by income level:
content(gas,"cons - mrio - avg int", fdtype, "lowincome2") = content(gas,"cons - direct - avg int", fdtype, "lowincome2") + sum((i,lowincome2), fd(fdtype,i,lowincome2)* co2int(gas,"MRIO cons indirect",i,"avg") )  / sum((i,lowincome2), fd(fdtype,i,lowincome2)) ;
content(gas,"cons - mrio - avg int", fdtype, "middleincome2") = content(gas,"cons - direct - avg int", fdtype, "middleincome2") + sum((i,middleincome2), fd(fdtype,i,middleincome2)* co2int(gas,"MRIO cons indirect",i,"avg") )  / sum((i,middleincome2), fd(fdtype,i,middleincome2)) ;
content(gas,"cons - mrio - avg int", fdtype, "lowincome") = content(gas,"cons - direct - avg int", fdtype, "lowincome") + sum((i,lowincome), fd(fdtype,i,lowincome)* co2int(gas,"MRIO cons indirect",i,"avg") )  / sum((i,lowincome), fd(fdtype,i,lowincome)) ;
content(gas,"cons - mrio - avg int", fdtype, "middleincome") = content(gas,"cons - direct - avg int", fdtype, "middleincome") + sum((i,middleincome), fd(fdtype,i,middleincome)* co2int(gas,"MRIO cons indirect",i,"avg") )  / sum((i,middleincome), fd(fdtype,i,middleincome)) ;
content(gas,"cons - mrio - avg int", fdtype, "highincome") = content(gas,"cons - direct - avg int", fdtype, "highincome") + sum((i,highincome), fd(fdtype,i,highincome)* co2int(gas,"MRIO cons indirect",i,"avg") )  / sum((i,highincome), fd(fdtype,i,highincome)) ;


* imported consumption 
content(gas,"imported cons - mrio", fdtype, r) = sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)*   trade_sh(i,s,r) * co2int(gas,"MRIO output",i,s))/sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)* trade_sh(i,s,r));
content(gas,"imported cons - mrio - avg int", fdtype, r) = sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)*   trade_sh(i,s,r) * co2int(gas,"MRIO output",i,"avg"))/sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)* trade_sh(i,s,r));

* as a function of total demand (income) 

content(gas,"imported cons - mrio (income)", fdtype, r) = sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)*   trade_sh(i,s,r) * co2int(gas,"MRIO output",i,s))/sum((i), fd(fdtype,i,r));
 
content(gas,"imported cons - mrio - avg int (income)", fdtype, r) = sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)*   trade_sh(i,s,r) * co2int(gas,"MRIO output",i,"avg"))/sum(i, fd(fdtype,i,r));


* share of imported co2 in consumed co2

content(gas,"imported cons share - mrio", fdtype, r)$sum(i, fd(fdtype,i,r)* co2int(gas,"MRIO cons total",i,r) ) = sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)*   trade_sh(i,s,r) * co2int(gas,"MRIO output",i,s))/  
 sum(i, fd(fdtype,i,r)* co2int(gas,"MRIO cons total",i,r) );

 
content(gas,"imported cons share - mrio - avg int", fdtype, r)$sum(i, fd(fdtype,i,r)* co2int(gas,"MRIO cons total",i,"avg") ) = sum((i,s)$(not sameas(r,s)), fd(fdtype,i,r)*   trade_sh(i,s,r) * co2int(gas,"MRIO output",i,"avg")) /  sum(i, fd(fdtype,i,r)* co2int(gas,"MRIO cons total",i,"avg") );

* share of indirect emissions in MRIO emissions
* note: want to put electricity in here

content(gas,"cons - share indirect", fdtype, r_) = 1-( content(gas,"cons - direct ele", fdtype, r_) / content(gas,"cons - mrio", fdtype, r_));

* at average int
content(gas,"cons - share indirect - avg int", fdtype, r_) =  1-( content(gas,"cons - direct ele - avg int", fdtype, r_) / content(gas,"cons - mrio - avg int", fdtype, r_));


* PRODUCTION 
* needs fitted production
* GETTING FITTED PRODUCTION FROM cf_compute_fittedprod_shares.gms

parameter fitted_production_absorption;

$gdxin estimates%SLASH%fitted_prod_shares_%spec%_%demandest%.gdx

$load fitted_production_absorption
fitted_production_absorption("production","data",i,r_) = fitted_production_absorption("production","true",i,r_);
 
* !! This does not include Final demand
content(gas,"prod (income)", fdtype, r_) = sum(i, co2int(gas,"output",i,r_) * fitted_production_absorption("production",fdtype,i,r_)) / sum((i), fd(fdtype,i,r_));

content(gas,"prod (income) - avg int", fdtype, r_) = sum(i, co2int(gas,"output",i,"avg") *fitted_production_absorption("production",fdtype,i,r_)) / sum((i), fd(fdtype,i,r_));


* as ratio of production
content(gas,"prod (prod)", fdtype, r_) = sum(i, co2int(gas,"output",i,r_) * fitted_production_absorption("production",fdtype,i,r_)) / sum(i, fitted_production_absorption("production",fdtype,i,r_));

content(gas,"prod (prod) - avg int", fdtype, r_) = sum(i, co2int(gas,"output",i,"avg") * fitted_production_absorption("production",fdtype,i,r_)) / sum(i, fitted_production_absorption("production",fdtype,i,r_));


* USE THIS
* computing holding technologies constant, j
* only variation is trade shares

content(gas,"prod (prod) - avg int ALT", fdtype, r_) = sum((s,i), co2int(gas,"MRIO output",i,"avg") * trade_sh(i,r_,s) * fd(fdtype,i,s)  ) / 
   sum((s,i),  trade_sh(i,r_,s) * fd(fdtype,i,s)  );
* denominator is production 

content(gas,"prod (prod) - avg int ALT", fdtype, "total") = sum((r_,s,i), co2int(gas,"MRIO output",i,"avg") * trade_sh(i,r_,s) * fd(fdtype,i,s)  ) ;
* check, totals are very similar: OK



content(gas,"prod (prod) - avg int ALT2", fdtype, r_) = sum((s,i), co2int(gas,"MRIO output",i,"avg") * trade_sh(i,r_,s) * fd(fdtype,i,s)  ) / 
  sum((i), fd(fdtype,i,r_));;

* total (= emissions occuring within country)
* as a function of income (~GDP)
* = prod + fd
content(gas,"total (income)", fdtype, r_) = (sum(i, co2int(gas,"output",i,r_) * fitted_production_absorption("production",fdtype,i,r_))  + sum(i_,co2coeff(i_,"fd",r_) * fd(fdtype,i_,r_))  )   / sum((i), fd(fdtype,i,r_));

content(gas,"total (income) - avg int", fdtype, r_) = (sum(i, co2int(gas,"output",i,"avg") * fitted_production_absorption("production",fdtype,i,r_))  + sum(i_,co2coeff(i_,"fd","avg") * fd(fdtype,i_,r_))  )   / sum((i), fd(fdtype,i,r_));


* Net exports ENX  
* = production emissions - embodied consumption emissions   (i.e. no final demand emissions)
content(gas,"net exports (income)", fdtype, r_) = (sum(i, co2int(gas,"output",i,r_) *fitted_production_absorption("production",fdtype,i,r_)) -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,r_) ))
		  / sum((i), fd(fdtype,i,r_)) ;


* take fitted production or not? doing so limits the differences between specifications

content(gas,"net exports (income) - avg int", fdtype, r_) = (sum(i, co2int(gas,"output",i,"AVG") * fitted_production_absorption("production",fdtype,i,r_)) -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") ))
		  / sum((i), fd(fdtype,i,r_)) ;

content(gas,"net exports (income)", fdtype, r_) = (sum(i, co2int(gas,"output",i,r_) *fitted_production_absorption("production",fdtype,i,r_)) -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,r_) ))
		  / sum((i), fd(fdtype,i,r_)) ;


* 
content(gas,"net exports (income) - avg int", fdtype, r_) = (sum(i, co2int(gas,"output",i,"AVG") * fitted_production_absorption("production",fdtype,i,r_)) -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") ))
		  / sum((i), fd(fdtype,i,r_)) ;

*  using alt production spec
* i.e. only diff between consumption and production is tradeshares..

content(gas,"net exports (income) - avg int ALT", fdtype, r_) = (  sum((s,i), co2int(gas,"MRIO output",i,"avg") * trade_sh(i,r_,s) * fd(fdtype,i,s)  )  -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") ))
		  / sum((i), fd(fdtype,i,r_)) ;

* not at intensity
content(gas,"net exports - avg int ALT", fdtype, r_) = (  sum((s,i), co2int(gas,"MRIO output",i,"avg") * trade_sh(i,r_,s) * fd(fdtype,i,s)  )  -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") ));

* use this to check sum
content(gas,"net exports - avg int ALT", fdtype, "total net exports") = sum(r_, content(gas,"net exports - avg int ALT", fdtype, r_) );


* as a function of consumption emissions -- BEST?
content(gas,"net exports (consumption) - avg int ALT", fdtype, r_) = (  sum((s,i), co2int(gas,"MRIO output",i,"avg") * trade_sh(i,r_,s) * fd(fdtype,i,s)  )  -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") ))
		  / sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") )  ;


* as a function of the emissions embodied in consumption

content(gas,"net exports (emissions)", fdtype, r_) = (sum(i, co2int(gas,"output",i,r_) *fitted_production_absorption("production",fdtype,i,r_)) -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,r_) ))
		  /  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,r_) ) ;

 
content(gas,"net exports (emissions) - avg int", fdtype, r_) = (sum(i, co2int(gas,"output",i,"AVG") * fitted_production_absorption("production",fdtype,i,r_)) -  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") ))
		  /  sum(i, fd(fdtype,i,r_)* co2int(gas,"MRIO cons indirect",i,"avg") ) ;



* ----------------------
* computing MODEL FIT


parameter	sse
		modelfit sum of square deviations to true values;
		
sse(gas,"avgcontent","cons - direct ele") =  sum(r_, sqr(content(gas,"cons - direct ele", "data", r_) - (sum(r_.local,content(gas,"cons - direct ele", "data", r_))/card(r_))));
sse(gas,"avgcontent","cons - mrio") =  sum(r_, sqr(content(gas,"cons - mrio", "data", r_) - (sum(r_.local,content(gas,"cons - mrio", "data", r_))/card(r_))));
sse(gas,"avgcontent","prod (prod)") =  sum(r_, sqr(content(gas,"prod (prod)", "data", r_) - (sum(r_.local,content(gas,"prod (prod)", "data", r_))/card(r_))));


modelfit(gas,"avgcontent","cons - direct ele",fdtype,"data")$sse(gas,"avgcontent","cons - direct ele") =
 1- sum(r_, sqr(content(gas,"cons - direct ele", fdtype, r_) - content(gas,"cons - direct ele", "data", r_))) / sse(gas,"avgcontent","cons - direct ele");


modelfit(gas,"avgcontent","cons - direct ele",fdtype,"avg")$sse(gas,"avgcontent","cons - direct ele") =
 1- sum(r_, sqr(content(gas,"cons - direct ele - avg int", fdtype, r_) - content(gas,"cons - direct ele", "data", r_)))/ sse(gas,"avgcontent","cons - direct ele");


modelfit(gas,"avgcontent","cons - mrio",fdtype,"data")$sse(gas,"avgcontent","cons - mrio") =
1- sum(r_, sqr(content(gas,"cons - mrio", fdtype, r_) - content(gas,"cons - mrio", "data", r_))) / sse(gas,"avgcontent","cons - mrio");


modelfit(gas,"avgcontent","cons - mrio",fdtype,"avg")$sse(gas,"avgcontent","cons - mrio") =
 1- sum(r_, sqr(content(gas,"cons - mrio - avg int", fdtype, r_) - content(gas,"cons - mrio", "data", r_)))/ sse(gas,"avgcontent","cons - mrio");


modelfit(gas,"avgcontent","prod (prod)",fdtype,"data")$sse(gas,"avgcontent","prod (prod)") =
1- sum(r_, sqr(content(gas,"prod (prod)", fdtype, r_) - content(gas,"prod (prod)", "data", r_))) / sse(gas,"avgcontent","prod (prod)");

modelfit(gas,"avgcontent","prod (prod)",fdtype,"avg")$sse(gas,"avgcontent","prod (prod)") =
 1- sum(r_, sqr(content(gas,"prod (prod) - avg int", fdtype, r_) - content(gas,"prod (prod)", "data", r_)))/ sse(gas,"avgcontent","prod (prod)");


* same thing weighted :


parameter weights;

weights(r_) = sum(i_, fd("data",i_,r_)) / sum((r_.local,i_), fd("data",i_,r_));


sse(gas,"avgcontent","cons - direct ele weighted") =  sum(r_, weights(r_) * sqr(content(gas,"cons - direct ele", "data", r_) - 
(sum(r_.local,sum(i_, fd("data",i_,r_)) * content(gas,"cons - direct ele", "data", r_))
/sum((r_.local,i_), fd("data",i_,r_)))));




sse(gas,"avgcontent","cons - mrio weighted") =  sum(r_, weights(r_) * sqr(content(gas,"cons - mrio", "data", r_) -
 (sum(r_.local,sum(i_, fd("data",i_,r_)) * content(gas,"cons - mrio", "data", r_))
/sum((r_.local,i_), fd("data",i_,r_)))));


sse(gas,"avgcontent","prod (prod) weighted") =  sum(r_, weights(r_) * sqr(content(gas,"prod (prod)", "data", r_) -
 (sum(r_.local,sum(i_, production(i_,r_)) * content(gas,"prod (prod)", "data", r_))
/sum((r_.local,i_), production(i_,r_)))));



modelfit(gas,"avgcontent weighted","cons - direct ele",fdtype,"data") =
 1- sum(r_, weights(r_) *sqr(content(gas,"cons - direct ele", fdtype, r_) - content(gas,"cons - direct ele", "data", r_))) / sse(gas,"avgcontent","cons - direct ele weighted");

modelfit(gas,"avgcontent weighted","cons - direct ele",fdtype,"avg") =
 1- sum(r_,weights(r_) * sqr(content(gas,"cons - direct ele - avg int", fdtype, r_) - content(gas,"cons - direct ele", "data", r_)))/ sse(gas,"avgcontent","cons - direct ele weighted");

modelfit(gas,"avgcontent weighted","cons - mrio",fdtype,"data") =
1- sum(r_, weights(r_) *sqr(content(gas,"cons - mrio", fdtype, r_) - content(gas,"cons - mrio", "data", r_))) / sse(gas,"avgcontent","cons - mrio weighted");

modelfit(gas,"avgcontent weighted","cons - mrio",fdtype,"avg") =
 1- sum(r_, weights(r_) *sqr(content(gas,"cons - mrio - avg int", fdtype, r_) - content(gas,"cons - mrio", "data", r_)))/ sse(gas,"avgcontent","cons - mrio weighted");

modelfit(gas,"avgcontent weighted","prod (prod)",fdtype,"data") =
1- sum(r_, weights(r_) *sqr(content(gas,"prod (prod)", fdtype, r_) - content(gas,"prod (prod)", "data", r_))) / sse(gas,"avgcontent","prod (prod) weighted");

modelfit(gas,"avgcontent weighted","prod (prod)",fdtype,"avg") =
 1- sum(r_, weights(r_) *sqr(content(gas,"prod (prod) - avg int", fdtype, r_) - content(gas,"prod (prod)", "data", r_)))/ sse(gas,"avgcontent","prod (prod) weighted");



* ---------------------------------------------------------------------
* (PARTIAL EQUILIBRIUM) APPROXIMATIONS:


* CO2/ ENERGY APPROXIMATIONS:
parameters  chg_fd_approx, chg_fd_qty_approx, chg_production_qty_approx,  chg_absorption_approx, chg_fprice_justfd, chg_fprice_approx, chg_skillprem_approx, chg_skillprem_justfd, log_Y_qty_approx;

** note: here "total demand elasticity" must be computed in datapreparation.gms !

*  p_c does not have a supply elasticity
supply_elast("p_c",r) = supply_elast("oil",r);

* for other sectors, set very high elasticity
supply_elast(i,r)$(not ff(i)) = 10000 ;


chg_fd_approx("approx",i,r) = 1 + (incelast("DIRECT",i,r) - 1) * (income_shock - 1);

chg_fd_approx("simulated",i,r) = chg_fd.L(i,r);


chg_fd_qty_approx("approx",i,r) = 1 + (incelast("DIRECT",i,r)) * (income_shock - 1);
chg_fd_qty_approx("simulated",i,r) = chg_fd.L(i,r) / chg_index.L(i,r);


chg_production_qty_approx("approx",i,r) = 1 + (incelast("total",i,r)) * (income_shock - 1);
chg_production_qty_approx("simulated",i,r) = chg_production.L(i,r) / chg_index.L(i,r);

* HOW WELL CAN WE APPROXIMATE CHANGES IN PRODUCTION?

log_Y_qty_approx("approx",i,r ) =  log(income_shock) * (incelast("total",i,r))  ;
log_Y_qty_approx("approx",pe,r ) =  log(income_shock) * (incelast("total",pe,r)) * (1 - 1 /(1+ supply_elast(pe,r))) ;

log_Y_qty_approx("simulated",i,r ) =  log( chg_production.L(i,r) / chg_Index.L(i,r)) ;


* INVERSION OF FINAL DEMAND VECTOR
* NOTE: this is now done USING TOTAL COEFFICIENTS and ARE COMPUTED IN DATAPREP.GMS FILE



* ---------------------------------------------------------------------
* SIMULATED CHANGES



parameter	chg_energy_co2,
		chg_totalfd_qty
		sh_co2_bilat,
		sh_co2_fd
		sh_co2_fd_ele
		sh_co2_fd_prod
		sh_co2_totalmrio
		sh_co2_region;

*		;


beta_secen(ene_all,i) = 1;
beta_secen(i, ene_all) = 0;


* co2

sh_co2_bilat("co2",ff,j,r) =  bmkco2(ff,j,r) / sum((ff.local,j.local), bmkco2(ff,j,r));  

sh_co2_fd("co2",ff,r) =  bmkco2(ff,"fd",r) / sum((ff.local), bmkco2(ff,"fd",r));  


* adding electricity
sh_co2_fd_ele("co2",e,r) =  co2coeff(e,"fd-ele",r) * fd("data",e,r) / sum(e.local, co2coeff(e,"fd-ele",r) * fd("data",e,r));

sh_co2_fd_prod("co2",r) =  sum(ff,bmkco2(ff,"fd",r)) / (sum((ff.local,j.local), bmkco2(ff,j,r)) + sum((ff.local), bmkco2(ff,"fd",r)));  

* secondary energy
sh_co2_bilat("secen",se,j,r) = energy("secondary",se,j,r) / sum((se.local,j.local), energy("secondary",se,j,r));  

sh_co2_fd("secen",se,r) = energy("secondary",se,"fd",r) / sum((se.local), energy("secondary",se,"fd",r));  

sh_co2_fd_prod("secen",r) =  sum(se,energy("secondary",se,"fd",r)) / (sum((se.local,j.local), energy("secondary",se,j,r)) + sum((se.local), energy("secondary",se,"fd",r)));  

* only for approx: 
sh_co2_totalmrio("co2", i,r) = fd("data",i,r)* co2int("co2","MRIO cons total",i,r) / sum(i.local, fd("data",i,r)* co2int("co2","MRIO cons total",i,r) ) ;


* regional shares for totals
* region's share in total emissions
sh_co2_region("co2 - total",r) =  (sum(ff,bmkco2(ff,"fd",r)) +  sum((ff.local,j.local), bmkco2(ff,j,r))) / sum(r.local, (sum(ff,bmkco2(ff,"fd",r)) +  sum((ff.local,j.local), bmkco2(ff,j,r))));
sh_co2_region("secen - total",r) =  (sum((se.local,j.local), energy("secondary",se,j,r)) + sum(se,energy("secondary",se,"fd",r))) / sum(r.LOCAL, (sum((se.local,j.local), energy("secondary",se,j,r)) + sum(se,energy("secondary",se,"fd",r))));

chg_energy_co2("logpci",r)  = log(pcexp(r));

* APPROXIMATIONS- CLOSED ECONOMY
* note, these imply no resprod

* note: here, ff includes p_c

* but p_c does not have a supply elasticity!

supply_elast("p_c",r) = supply_elast("oil",r);

* P.E. APPROXIMATIONS USED FOR NOW
* ignoring the effect of the supply elasticity
chg_energy_co2("co2 - direct - approx simple",r) = sum(ff, sh_co2_FD("co2",ff,r) * ( incelast("direct",ff,r)));
chg_energy_co2("co2 - direct ele - approx simple",r) = sum(se, sh_co2_fd_ele("co2",se,r) * ( incelast("direct",se,r)));
chg_energy_co2("co2 - total - approx simple",r) = sum(i, sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,r)));
chg_energy_co2("co2 - total - approx full",r) = sum(i, sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,r))) ;


* with the effect of supply elasticity:

* ignore these:
chg_energy_co2("co2 - prod - approx",r) = exp( ( sum(j, sum(ff, sh_co2_bilat("co2",ff,j,r)) * incelast("total",j,r)) - sum(ff, sum(j, sh_co2_bilat("co2",ff,j,r)) * incelast("total",ff,r) / (1 + supply_elast(ff,r) ))) * log(income_shock));

chg_energy_co2("co2 - direct - approx",r) = exp(  sum(ff, sh_co2_FD("co2",ff,r) * ( incelast("direct",ff,r) - incelast("total",ff,r) / (1 + supply_elast(ff,r) ))) * log(income_shock));

chg_energy_co2("co2 - total - approx",r) = sh_co2_fd_prod("co2",r) * chg_energy_co2("co2 - direct - approx",r) + (1-sh_co2_fd_prod("co2",r)) * chg_energy_co2("co2 - prod - approx",r) ;


*  what to do with electricity?
*using supply elasticity  of fossil fuels:
supply_elast("ely",r) = supply_elast("oil",r);

* with no trade
chg_energy_co2("co2 - direct ele - approx full no trade",r) = sum(se, sh_co2_fd_ele("co2",se,r) * ( incelast("direct",se,r))) - 
   sum(se, sh_co2_fd_ele("co2",se,r) * incelast("total",se,r) / (1 + supply_elast(se,r) ));


* with trade 
chg_energy_co2("co2 - direct ele - approx full",r) = sum(se, sh_co2_fd_ele("co2",se,r) * ( incelast("direct",se,r))) - 
   sum(se, sh_co2_fd_ele("co2",se,r) *
   sum(s, trade_sh(se,s,r) * incelast("total",se,s) / (1 + supply_elast(se,s) )));


* decomposition

set decompset / "co2 - total" "co2 - total USA share", "co2 - total USA incel", "co2 - total avg coeffs", "co2 - total - sum zero", "co2 - total USA share - sum zero","co2 - total USA incel - sum zero",
"co2 - total avg coeffs - sum zero" ,"co2 - direct ele" ,"co2 - direct ele - sum zero", "co2 - total - simul fd - sum zero", "FOR TABLE - co2 - total - sum zero"/;

parameter decomposition;
decomposition("co2 - total",i,r) =  sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,r)   );

decomposition("co2 - total USA share",i,r) =  sh_co2_totalmrio("co2", i,"usa") * ( incelast("direct",i,r)   );
decomposition("co2 - total USA incel",i,r) =  sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,"usa")   );

* using average co2 coeff
decomposition("co2 - total avg coeffs",i,r) = fd("data",i,r)* co2int("co2","MRIO cons total",i,"avg") / sum(i.local, fd("data",i,r)* co2int("co2","MRIO cons total",i,"avg"))  * ( incelast("direct",i,r)   );

decomposition("co2 - total - sum zero",i,r) =  sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,r) -1  );


decomposition("co2 - total - sum zero",i,"world") = sum((r), sh_co2_region("co2 - total",r) * decomposition("co2 - total - sum zero",i,r));


* using simulated change:
decomposition("co2 - total - simul fd - sum zero",i,r) =  sh_co2_totalmrio("co2", i,r) * ( chg_fd.L(i,r) / chg_Index.L(i,r) -1   ) / (income_shock - 1)    ;


decomposition("co2 - total USA share - sum zero",i,r) =  sh_co2_totalmrio("co2", i,"usa") * ( incelast("direct",i,r) -1  );
decomposition("co2 - total USA incel - sum zero",i,r) =  sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,"usa")-1   );

* using average co2 coeff
decomposition("co2 - total avg coeffs - sum zero",i,r) = fd("data",i,r)* co2int("co2","MRIO cons total",i,"avg") / sum(i.local, fd("data",i,r)* co2int("co2","MRIO cons total",i,"avg"))  * ( incelast("direct",i,r)  -1 );

decomposition("co2 - direct ele",e,r) =  sh_co2_fd_ele("co2",e,r) * ( incelast("direct",e,r));

decomposition("co2 - direct ele - sum zero",e,r) =  sh_co2_fd_ele("co2",e,r) * ( incelast("direct",e,r));


decomposition(decompset,"logpci",r) =  log(pcexp(r));
decomposition(decompset,"total",r) =  sum(i, decomposition(decompset,i,r));


* FOR PAPER - DECOMPOSITION TABLE

decomposition("FOR TABLE - co2 - total - sum zero","total",r) = sum(i, sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,r) -1  ));

decomposition("FOR TABLE - co2 - total - sum zero",se,r) =  sh_co2_totalmrio("co2", se,r) * ( incelast("direct",se,r) -1  );
decomposition("FOR TABLE - co2 - total - sum zero","energy goods",r) = sum(se, sh_co2_totalmrio("co2", se,r) * ( incelast("direct",se,r) -1  ));

* aggregate to large sectors 
decomposition("FOR TABLE - co2 - total - sum zero",ii_,r) =  sum(i$mapi(i,ii_), sh_co2_totalmrio("co2", i,r) * ( incelast("direct",i,r) -1  ));


decomposition("FOR TABLE - co2 - total - sum zero","non-energy goods",r) = decomposition("FOR TABLE - co2 - total - sum zero","total",r) - 
decomposition("FOR TABLE - co2 - total - sum zero","energy goods",r);

* aggregate across countries
* make a set of aggregation
set aggregset / "total", "energy goods", "non-energy goods" /;
aggregset(se) = yes;
aggregset(ii_) = yes;

* using new definition:

decomposition("FOR TABLE - co2 - total - sum zero",aggregset,"lowincome") = sum(lowincome, sh_co2_region("co2 - total",lowincome) * decomposition("FOR TABLE - co2 - total - sum zero",aggregset,lowincome)) /
			sum(lowincome,  sh_co2_region("co2 - total",lowincome)) ;


decomposition("FOR TABLE - co2 - total - sum zero",aggregset,"lowincome2") = sum(lowincome2, sh_co2_region("co2 - total",lowincome2) * decomposition("FOR TABLE - co2 - total - sum zero",aggregset,lowincome2)) /
			sum(lowincome2,  sh_co2_region("co2 - total",lowincome2)) ;


decomposition("FOR TABLE - co2 - total - sum zero",aggregset,"middleincome") = sum( middleincome, sh_co2_region("co2 - total", middleincome) * decomposition("FOR TABLE - co2 - total - sum zero",aggregset, middleincome)) /
			sum( middleincome,  sh_co2_region("co2 - total", middleincome)) ;

decomposition("FOR TABLE - co2 - total - sum zero",aggregset,"middleincome2") = sum( middleincome2, sh_co2_region("co2 - total", middleincome2) * decomposition("FOR TABLE - co2 - total - sum zero",aggregset, middleincome2)) /
			sum( middleincome2,  sh_co2_region("co2 - total", middleincome2)) ;

decomposition("FOR TABLE - co2 - total - sum zero",aggregset,"highincome") = sum( highincome, sh_co2_region("co2 - total", highincome) * decomposition("FOR TABLE - co2 - total - sum zero",aggregset, highincome)) /
			sum( highincome,  sh_co2_region("co2 - total", highincome)) ;

decomposition("FOR TABLE - co2 - total - sum zero",aggregset,"world") = sum((r), sh_co2_region("co2 - total",r) * decomposition("FOR TABLE - co2 - total - sum zero",aggregset,r));


* ----------------------
* EXACT (SIMULATED) CHANGES

* NOTE: this computes direct, direct-ele and production emissions ("total"), but the change in MRIO emissions is harder to estimate and is done by comparing consumption contents
* before and after shock


chg_energy_co2("co2 - prod",r) = sum((i,j), sh_co2_bilat("co2",i,j,r) * chg_production.L(j,r)   /  chg_Index.L(i,r));

* in final demand
chg_energy_co2("co2 - direct",r) = sum((i), sh_co2_fd("co2",i,r) * chg_fd.L(i,r)   /  chg_Index.L(i,r));

* direct ELE is computed further down -- SINCE IT NEEDS MRIO INTENSITY

chg_energy_co2("co2 - total",r) = sh_co2_fd_prod("co2",r) * chg_energy_co2("co2 - direct",r) + (1-sh_co2_fd_prod("co2",r)) * chg_energy_co2("co2 - prod",r) ;

chg_energy_co2("co2 - total","world") = sum(r, sh_co2_region("co2 - total",r) * chg_energy_co2("co2 - total",r));


chg_energy_co2("secen - prod",r) = sum((i,j), sh_co2_bilat("secen",i,j,r) * chg_production.L(j,r)   /  chg_Index.L(i,r));

chg_energy_co2("secen - direct",r) = sum((i), sh_co2_fd("secen",i,r) * chg_fd.L(i,r)   /  chg_Index.L(i,r));

chg_energy_co2("secen - total",r) = sh_co2_fd_prod("secen",r) * chg_energy_co2("secen - direct",r) + (1-sh_co2_fd_prod("secen",r)) * chg_energy_co2("secen - prod",r) ;

chg_energy_co2("secen - total","world") = sum(r, sh_co2_region("secen - total",r) * chg_energy_co2("secen - total",r));


* intensities: 

content(gas,"income", "data", r) = sum(i,     fd("data",i,r));
content(gas,"income - counterfact", "data", r) = sum(i,  chg_fd.L(i,r)  /  chg_Index.L(i,r)  *   fd("data",i,r));

* here CHANGE as function of total change in final demand (qty terms)
* This correspongd to change in real income

* equiv to difference in income

chg_totalfd_qty(r) = sum((i), dem_sh(i,r)  * chg_fd.L(i,r)   /  chg_Index.L(i,r));

chg_totalfd_qty("world") = sum((r,i),fd("data",i,r) * chg_totalfd_qty(r))  / sum((i.local,r), fd("data",i,r));


* change in production
chg_production_qty(r) = sum(i, chg_production.L(i,r)   /  chg_Index.L(i,r)  * 
  (fitted_production_absorption("production","non-homoth",i,r) / sum(i.local, fitted_production_absorption("production","non-homoth",i,r)   )));


* this is production
chg_energy_co2("co2 - int - total",r) = chg_energy_co2("co2 - total",r) / chg_totalfd_qty(r);
chg_energy_co2("secen - int - total",r) = chg_energy_co2("secen - total",r) / chg_totalfd_qty(r);

* as elasticity 

chg_energy_co2("co2 - elast - direct",r) = (chg_energy_co2("co2 - direct",r) -1) / (chg_totalfd_qty(r) -1);

* production only -- to compare to MRIO indirect
chg_energy_co2("co2 - elast - prod",r) = (chg_energy_co2("co2 - prod",r) -1) / (chg_totalfd_qty(r) -1);

chg_energy_co2("co2 - elast - total",r) = (chg_energy_co2("co2 - total",r) -1) / (chg_totalfd_qty(r) -1);

* USING CHG IN PRODUCTION AS DENOMINATOR
chg_energy_co2("co2 - elast - total - denom prod",r) = (chg_energy_co2("co2 - total",r) -1) / (chg_production_qty(r)  -1);

chg_energy_co2("secen - elast - total",r) = (chg_energy_co2("secen - total",r) -1) / (chg_totalfd_qty(r) -1);

chg_energy_co2("co2 - elast - total","world") = (chg_energy_co2("co2 - total","world") -1) / (chg_totalfd_qty("world") -1);
chg_energy_co2("secen - elast - total","world") = (chg_energy_co2("secen - total","world") -1) / (chg_totalfd_qty("world") -1);


* EMBODIED EMISSIONS:
* requires re-computing the MRIO coefficients


* after simultaion import shares:
parameter bitc_shock;
bitc_shock(r,s) = 1;

parameter trade_sh_counterfactual(i,*,s), production_counterfactual;

*
trade_sh_counterfactual(i,r,s) = trade_sh(i,r,s) * chg_supplier.L(i,r) * bitc_shock(r,s)**(-theta(i)) * chg_index.L(i,s)**theta(i);

* just to check -- note, trade_sh don't add up to one exactly
trade_sh_counterfactual(i,"total",s) = sum(r, trade_sh_counterfactual(i,r,s))  ;

trade_sh_counterfactual(i,r,s)$trade_sh_counterfactual(i,"total",s) = trade_sh_counterfactual(i,r,s) / trade_sh_counterfactual(i,"total",s);


* !!! iter MUST be the same as in other file
set	iter /0*10/;

parameters	ay(*,i,r)		MRIO carbon content of output
		ay2(i,r)	to compute SECONDARY energy content
*		af(g,r)		Direct carbon content (fuel only)
*		ad(g,r)		Direct carbon content (fuel + electricity)
		am(*,i,r)		Carbon content of imports
		am2
		at(j)		Carbon content of trade,
		at2
		atf(j)		Carbon content of trade (fuel only),
		atd(j)		Carbon content of trade (fuel + electricity)
		co2e(*,*,r)	Total carbon emissions by sector and region;



co2e("co2",j,r) = sum(i,bmkco2(i,j,r) *  chg_production.L(j,r)   /  chg_Index.L(i,r));

production_counterfactual(i,r)  = production(i,r) *  chg_production.L(i,r)   /  chg_Index.L(i,r);

* do them all
*set gas_ /co2, nh4, n20, fgas, ghg/;

* only co2 and total GHG
set gas_ /co2, ghg/;



* -- !! TO MAKE SURE WE ARE COMPARING THE COMPARABLE, RUN THIS PLACEBO TEST:
* as placebo, need to run without change:
*co2e(j,r) = sum(i,bmkco2(i,j,r));
*trade_sh_counterfactual(i,r,s) = trade_sh(i,r,s);
*production_counterfactual(i,r) = production(i,r);


* OK!
*chg_fd.L(i,r) = 1;

ay(gas_,j,r)$production(j,r) = co2e(gas_,j,r) / production_counterfactual(j,r);


ay2(j,r)$production(j,r) = sum(i,energy("secondary",i,j,r)) / production_counterfactual(j,r);

parameter	dev		Maximum deviation in embodied carbon
		aynew		New estimate of embodied carbon,
		ay2new
*		ab(i,r,s,*,metric)	Carbon content of bilateral trade flow (gross and net);
;

loop(iter,
	aynew(gas_,i,r)$production(i,r) =  co2e(gas_,i,r) / production_counterfactual(i,r) + sum((j,s), trade_sh_counterfactual(j,s,r) *  intersh(j,i,r) * ay(gas_,j,s))	;

	dev(gas_,iter,r) = smax(i, abs(ay(gas_,i,r)-aynew(gas_,i,r)));
	ay(gas_,i,r) = aynew(gas_,i,r);


* now for energy

	ay2new(i,r)$production(i,r) =  sum(j,energy("secondary",j,i,r)) / production_counterfactual(i,r)+ sum((j,s), trade_sh_counterfactual(j,s,r) *  intersh(j,i,r) * ay2(j,s))	;
	ay2(i,r) = ay2new(i,r);

);

* import intensity
am(gas_,i,r)$sum(s$(not sameas(r,s)),  trade_sh(i,s,r)) =  sum(s$(not sameas(r,s)),  trade_sh_counterfactual(i,s,r) *ay(gas_,i,s))/sum(s$(not sameas(r,s)),  trade_sh(i,s,r));

am2(i,r)$sum(s$(not sameas(r,s)),  trade_sh(i,s,r)) =  sum(s$(not sameas(r,s)),  trade_sh_counterfactual(i,s,r) *ay2(i,s))/sum(s$(not sameas(r,s)),  trade_sh(i,s,r));


co2int(gas_,"MRIO output - counterfact",i,r) = ay(gas_,i,r);

co2int(gas_,"MRIO cons indirect - counterfact",i,r) = sum(s,  trade_sh_counterfactual(i,s,r) *ay(gas_,i,s));

co2int(gas_,"MRIO imports - counterfact",i,r) = am(gas_,i,r);

* computing the counterfactual change in embodied emissions
* this is total: MRIO + direct

content("co2","cons - direct - counterfact", "data", r) =   sum(i,co2coeff(i,"fd",r) * chg_fd.L(i,r)   /  chg_Index.L(i,r) * fd("data",i,r)) / content("co2","income - counterfact", "data", r);

content("co2","cons - direct ele - counterfact", "data", r) =   sum(i,co2coeff(i,"fd-ele",r) * chg_fd.L(i,r)   /  chg_Index.L(i,r) * fd("data",i,r)) / content("co2","income - counterfact", "data", r);

* should be in qty.. correct?
content("co2","cons - mrio - counterfact", "data", r) = 
	content("co2","cons - direct - counterfact", "data", r) +
	 sum(i, fd("data",i,r) * chg_fd.L(i,r) / chg_Index.L(i,r)  * co2int("co2","MRIO cons indirect - counterfact",i,r) )
		  / content("co2","income - counterfact", "data", r) ;

* indirect only (for net exports)
content("co2","cons - mrio indirect - counterfact", "data", r) = 
	 sum(i, fd("data",i,r) * chg_fd.L(i,r) / chg_Index.L(i,r)  * co2int("co2","MRIO cons indirect - counterfact",i,r) )
		  / content("co2","income - counterfact", "data", r) ;

* this is ECI_hat

*PROBABLY better to MAKE A NEW PARAMETER FOR THINGS THAT ARE NOT DEFINED AS INTENSITIES?

* multiply by income to get total change
chg_energy_co2("co2 - cons - mrio indirect",r) =  (content("co2","cons - mrio indirect - counterfact", "data", r) * content("co2","income - counterfact", "data", r)
	 - content("co2","cons - mrio indirect", "data", r)  * content("co2","income", "data", r)    ) /
	( content("co2","cons - mrio indirect", "data", r) * content("co2","income", "data", r)) + 1;

* elasticity
chg_energy_co2("co2 - elast - cons - mrio indirect",r) = (content("co2","cons - mrio indirect - counterfact", "data", r) * content("co2","income - counterfact", "data", r)
	 / content("co2","cons - mrio indirect", "data", r)  / content("co2","income", "data", r)   -1) 
			/ (chg_totalfd_qty(r)  -1);



* compute the elasticity as the pct difference between the original and counterfactual embodied co2 in consumption

* note, here computing as ratio of intensities, but what I really want is ratio of actual contents (hence multiplying by income change)

chg_energy_co2("co2 - elast - cons - direct 2",r) = (content("co2","cons - direct - counterfact", "data", r) * content("co2","income - counterfact", "data", r) / content("co2","cons - direct", "data", r) 
	/ content("co2","income", "data", r)   -1) 
			/ (chg_totalfd_qty(r)  -1);




chg_energy_co2("co2 - elast - cons - direct ele",r) = (content("co2","cons - direct ele - counterfact", "data", r) * content("co2","income - counterfact", "data", r) / content("co2","cons - direct ele", "data", r) 
	/ content("co2","income", "data", r)   -1) 
			/ (chg_totalfd_qty(r)  -1);


chg_energy_co2("co2 - elast - cons - direct ele","world") = (sum(r, content("co2","cons - direct ele - counterfact", "data", r) * content("co2","income - counterfact", "data", r)  )
	 / sum(r,content("co2","cons - direct ele", "data", r) * content("co2","income", "data", r)) -1) 
			/ (chg_totalfd_qty("world") -1);




chg_energy_co2("co2 - elast - cons mrio",r) = (content("co2","cons - mrio - counterfact", "data", r) * content("co2","income - counterfact", "data", r) / content("co2","cons - mrio", "data", r) 
	/ content("co2","income", "data", r)   -1) 
			/ (chg_totalfd_qty(r)  -1);


chg_energy_co2("co2 - elast - cons mrio","world") = (sum(r, content("co2","cons - mrio - counterfact", "data", r) * content("co2","income - counterfact", "data", r)  )
	 / sum(r,content("co2","cons - mrio", "data", r) * content("co2","income", "data", r)) -1) 
			/ (chg_totalfd_qty("world")  -1);


* note: the 'content' parameters are actually intensities.. need to multiply by income


chg_energy_co2("co2 - change - cons - direct ele",r) =  (content("co2","cons - direct ele - counterfact", "data", r) * content("co2","income - counterfact", "data", r) 
	 - content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) ) 
	/ (content("co2","cons - direct ele", "data", r ) * content("co2","income", "data", r));

chg_energy_co2("co2 - change - cons - mrio",r) = (content("co2","cons - mrio - counterfact", "data", r)* content("co2","income - counterfact", "data", r)  - content("co2","cons - mrio", "data", r) * content("co2","income", "data", r) ) 
	/ (content("co2","cons - mrio", "data", r )*content("co2","income", "data", r));


* compare to actual income growth:
chg_energy_co2("co2 - change - income",r) = ( content("co2","income - counterfact", "data", r)  -  content("co2","income", "data", r) ) 
	/ (content("co2","income", "data", r));

chg_energy_co2("co2 - change - income","world") = sum(r,(content("co2","income - counterfact", "data", r)  -  content("co2","income", "data", r) )) 
	/ sum(r,(content("co2","income", "data", r)));

* .. is pretty close. OK

chg_energy_co2("co2 - change - cons - direct ele","world") = sum(r, (content("co2","cons - direct ele - counterfact", "data", r) * content("co2","income - counterfact", "data", r) -
    content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)  )) 
	/ sum(r,content("co2","cons - direct ele", "data", r )*content("co2","income", "data", r)   );


chg_energy_co2("co2 - change - cons - mrio","world") = sum(r, (content("co2","cons - mrio - counterfact", "data", r) * content("co2","income - counterfact", "data", r) - content("co2","cons - mrio", "data", r)*content("co2","income", "data", r)  )) 
	/ sum(r,content("co2","cons - mrio", "data", r ) *content("co2","income", "data", r) );


* as elasticities:
chg_energy_co2("co2 - as elast - cons - direct ele",r)= chg_energy_co2("co2 - change - cons - direct ele",r) / chg_energy_co2("co2 - change - income",r);
chg_energy_co2("co2 - as elast - cons - mrio",r)= chg_energy_co2("co2 - change - cons - mrio",r) / chg_energy_co2("co2 - change - income",r);

chg_energy_co2("co2 - as elast - cons - direct ele","world")= chg_energy_co2("co2 - change - cons - direct ele","world") / chg_energy_co2("co2 - change - income","world");
chg_energy_co2("co2 - as elast - cons - mrio","world")= chg_energy_co2("co2 - change - cons - mrio","world") / chg_energy_co2("co2 - change - income","world");


*--------------------


* here, do the same for low, mid and high inc countries

chg_totalfd_qty("lowincome2") = sum((lowincome2,i),fd("data",i,lowincome2) * chg_totalfd_qty(lowincome2))  / sum((i.local,lowincome2), fd("data",i,lowincome2));

chg_energy_co2("co2 - elast - cons - direct ele","lowincome2") = (sum(lowincome2, content("co2","cons - direct ele - counterfact", "data", lowincome2) * content("co2","income - counterfact", "data", lowincome2)  )
	 / sum(lowincome2,content("co2","cons - direct ele", "data", lowincome2) * content("co2","income", "data", lowincome2)) -1) 
			/ (chg_totalfd_qty("lowincome2") -1);

chg_energy_co2("co2 - elast - cons mrio","lowincome2") = (sum(lowincome2, content("co2","cons - mrio - counterfact", "data", lowincome2) * content("co2","income - counterfact", "data", lowincome2)  )
	 / sum(lowincome2,content("co2","cons - mrio", "data", lowincome2) * content("co2","income", "data", lowincome2)) -1) 
			/ (chg_totalfd_qty("lowincome2")  -1);


chg_totalfd_qty("middleincome2") = sum((middleincome2,i),fd("data",i,middleincome2) * chg_totalfd_qty(middleincome2))  / sum((i.local,middleincome2), fd("data",i,middleincome2));

chg_energy_co2("co2 - elast - cons - direct ele","middleincome2") = (sum(middleincome2, content("co2","cons - direct ele - counterfact", "data", middleincome2) * content("co2","income - counterfact", "data", middleincome2)  )
	 / sum(middleincome2,content("co2","cons - direct ele", "data", middleincome2) * content("co2","income", "data", middleincome2)) -1) 
			/ (chg_totalfd_qty("middleincome2") -1);

chg_energy_co2("co2 - elast - cons mrio","middleincome2") = (sum(middleincome2, content("co2","cons - mrio - counterfact", "data", middleincome2) * content("co2","income - counterfact", "data", middleincome2)  )
	 / sum(middleincome2,content("co2","cons - mrio", "data", middleincome2) * content("co2","income", "data", middleincome2)) -1) 
			/ (chg_totalfd_qty("middleincome2")  -1);



chg_totalfd_qty("lowincome") = sum((lowincome,i),fd("data",i,lowincome) * chg_totalfd_qty(lowincome))  / sum((i.local,lowincome), fd("data",i,lowincome));

chg_energy_co2("co2 - elast - cons - direct ele","lowincome") = (sum(lowincome, content("co2","cons - direct ele - counterfact", "data", lowincome) * content("co2","income - counterfact", "data", lowincome)  )
	 / sum(lowincome,content("co2","cons - direct ele", "data", lowincome) * content("co2","income", "data", lowincome)) -1) 
			/ (chg_totalfd_qty("lowincome") -1);

chg_energy_co2("co2 - elast - cons mrio","lowincome") = (sum(lowincome, content("co2","cons - mrio - counterfact", "data", lowincome) * content("co2","income - counterfact", "data", lowincome)  )
	 / sum(lowincome,content("co2","cons - mrio", "data", lowincome) * content("co2","income", "data", lowincome)) -1) 
			/ (chg_totalfd_qty("lowincome")  -1);


chg_totalfd_qty("middleincome") = sum((middleincome,i),fd("data",i,middleincome) * chg_totalfd_qty(middleincome))  / sum((i.local,middleincome), fd("data",i,middleincome));

chg_energy_co2("co2 - elast - cons - direct ele","middleincome") = (sum(middleincome, content("co2","cons - direct ele - counterfact", "data", middleincome) * content("co2","income - counterfact", "data", middleincome)  )
	 / sum(middleincome,content("co2","cons - direct ele", "data", middleincome) * content("co2","income", "data", middleincome)) -1) 
			/ (chg_totalfd_qty("middleincome") -1);

chg_energy_co2("co2 - elast - cons mrio","middleincome") = (sum(middleincome, content("co2","cons - mrio - counterfact", "data", middleincome) * content("co2","income - counterfact", "data", middleincome)  )
	 / sum(middleincome,content("co2","cons - mrio", "data", middleincome) * content("co2","income", "data", middleincome)) -1) 
			/ (chg_totalfd_qty("middleincome")  -1);



chg_totalfd_qty("highincome") = sum((highincome,i),fd("data",i,highincome) * chg_totalfd_qty(highincome))  / sum((i.local,highincome), fd("data",i,highincome));

chg_energy_co2("co2 - elast - cons - direct ele","highincome") = (sum(highincome, content("co2","cons - direct ele - counterfact", "data", highincome) * content("co2","income - counterfact", "data", highincome)  )
	 / sum(highincome,content("co2","cons - direct ele", "data", highincome) * content("co2","income", "data", highincome)) -1) 
			/ (chg_totalfd_qty("highincome") -1);

chg_energy_co2("co2 - elast - cons mrio","highincome") = (sum(highincome, content("co2","cons - mrio - counterfact", "data", highincome) * content("co2","income - counterfact", "data", highincome)  )
	 / sum(highincome,content("co2","cons - mrio", "data", highincome) * content("co2","income", "data", highincome)) -1) 
	/ (chg_totalfd_qty("highincome")  -1);


* alt denom -- using actual shock size. 
*-- does not help --
chg_energy_co2("co2 - elast - cons mrio - alt",r) = (content("co2","cons - mrio - counterfact", "data", r) * content("co2","income - counterfact", "data", r) / content("co2","cons - mrio", "data", r) 
	/ content("co2","income", "data", r)   -1) 
			/ (1.01  -1);


chg_energy_co2("co2 - elast - cons mrio - alt","world") = (sum(r, content("co2","cons - mrio - counterfact", "data", r) * content("co2","income - counterfact", "data", r)  )
	 / sum(r,content("co2","cons - mrio", "data", r) * content("co2","income", "data", r)) -1) 
			/ (1.01  -1);

* NET EXPORTS: if
* better to define as ratio. But then what is interpretation of change ?
* net exports go up by 

* as income grows by 1%, the ratio of 
* need to compute indirect only ?


* as a ratio
chg_energy_co2("co2 - elast - net exports ratio",r)  = (chg_energy_co2("co2 - prod",r)  / chg_energy_co2("co2 - cons - mrio indirect",r) -1) / (chg_totalfd_qty(r)  -1);


*---------------------------------------------

* FOR REVISION, WHEN USING OBSERVED GDP GROWTH

* if not for simulation, simply an ex-post evaluation..

* import p.c. GDP numbers from 
* GAMS\data\pcgdpe_1990_2005_2014.dta

* PAST GROWTH

table pcgdp 
	1990	2005	2014
AGO	3176	3736	7860
ALB	3796	6104	10581
ARG	5931	14264	20144
ARM	5669	4979	8576
ATG	11920	17374	20817
AUS	26759	40212	42663
AUT	25245	38059	47285
AZE	6967	4652	15844
BDI	749	551	771
BEL	25817	36109	43203
BEN	1181	1734	1909
BFA	890	1229	1550
BGD	1379	1492	2862
BGR	11068	10638	17359
BHR	15485	38968	41325
BHS	20498	28082	23166
BIH	1304	6483	9976
BLR	12591	10304	20177
BLZ	5844	7407	8377
BMU	28811	49854	57193
BOL	2115	3735	5988
BRA	5678	8732	14742
BRB	16209	20638	14126
BRN	52401	71443	67558
BTN	2382	5247	6751
BWA	6214	12836	15891
CAF	1046	843	594
CAN	30293	40770	41936
CHE	38881	43874	57855
CHL	7638	13623	21404
CHN	2363	6520	12240
CIV	2263	2046	3336
CMR	2348	2305	2673
COD	1216	554	1211
COG	1858	3233	4354
COL	6860	7643	12516
COM	1838	1465	1460
CPV	2085	4181	6212
CRI	7221	10393	14099
CYP	21036	30488	28437
CZE	21267	24318	31567
DEU	25465	35730	45652
DJI	2496	2516	3162
DMA	7244	8132	10134
DNK	25970	36768	44522
DOM	4701	8657	12445
DZA	7882	11520	12650
ECU	4993	7317	10871
EGY	2031	5308	9893
ESP	17693	30536	33568
EST	10860	17808	28247
ETH	848	630	1308
FIN	24976	35222	40027
FJI	5649	6144	7858
FRA	24742	33357	39029
GAB	8436	17532	13990
GBR	24327	38251	39960
GEO	10201	4349	9332
GHA	1886	2450	3531
GIN	2088	1074	1427
GMB	1696	1721	1539
GNB	1253	1289	1252
GNQ	510	18729	39411
GRC	16336	27387	25888
GRD	4558	9896	11104
GTM	3610	5218	6835
HKG	24870	44499	51400
HND	2859	3642	4397
HRV	13836	16986	21519
HUN	12149	18818	25586
IDN	3217	4321	9576
IND	1310	2748	5168
IRL	18087	43096	48283
IRN	2976	14166	15409
IRQ	4395	4359	12033
ISL	30515	42621	42654
ISR	20942	29263	33024
ITA	25384	32683	35531
JAM	4449	6717	7402
JOR	3551	5342	10397
JPN	26714	34221	35107
KAZ	10812	10249	23340
KEN	2032	1885	2757
KGZ	6978	2052	3350
KHM	1051	1900	2983
KNA	9664	17073	23120
KOR	12347	28597	34701
KWT	20055	65456	63486
LAO	1297	2145	5491
LBN	3841	12887	13852
LBR	522	512	834
LCA	6985	9192	9990
LKA	2825	4831	10262
LSO	1116	1868	2390
LTU	12066	15900	28048
LUX	44423	74931	94224
LVA	15274	14941	23487
MAC	23679	45532	126216
MAR	4008	4400	7063
MDA	5042	2558	4805
MDG	990	1022	1227
MDV	4203	9373	14244
MEX	10208	13318	15745
MKD	6511	8435	13032
MLI	701	1154	1427
MLT	13611	24648	31432
MMR	790	1854	5299
MNE	9686	8966	14494
MNG	3141	4919	11321
MOZ	509	907	1118
MRT	1864	2212	3357
MUS	9140	12489	17777
MWI	947	1007	949
MYS	8266	16179	22949
NAM	4919	6453	10783
NER	823	706	843
NGA	844	3436	5485
NLD	25657	40215	46847
NOR	26405	51020	63479
NPL	1000	1291	2155
NZL	20642	28759	34499
OMN	10605	27373	38075
PAK	2443	3139	4641
PAN	6052	11815	19437
PER	3395	6435	10915
PHL	3452	4062	6624
POL	7771	14948	25028
PRT	14426	24157	28306
PRY	3686	4825	8258
QAT	26221	106209	142044
ROU	6479	10106	20636
RUS	18323	13757	23918
RWA	997	990	1557
SAU	18588	27828	47441
SDN	1600	2341	3764
SEN	1992	1936	2228
SGP	20767	57468	71568
SLE	1541	1059	1411
SLV	1022	3487	7834
SRB	7256	9013	13376
STP	1842	2131	3185
SUR	6502	9912	15421
SVK	16171	17943	28421
SVN	18940	26447	30250
SWE	27115	36705	44175
SWZ	5697	6979	7999
SYR	1123	5098	4187
TCD	1330	1853	2006
TGO	1333	1084	1377
THA	5015	9420	13843
TJK	5892	2057	2745
TKM	9677	9510	20678
TTO	11177	22138	31052
TUN	5208	8485	10286
TUR	9887	12355	19099
TWN	18852	36066	43999
TZA	1116	1477	2189
UGA	801	1303	1824
UKR	10218	7036	10324
URY	8936	10299	20244
USA	36620	50783	51983
UZB	4215	4613	8165
VCT	5194	8331	9542
VEN	8438	10178	14007
VNM	1272	2984	5310
YEM	662	3156	3342
ZAF	8211	10245	12044
ZMB	1354	1585	3672
ZWE	4148	1317	1867;

display pcgdp;


* DOING BOTH 1990 AND 2005 AS STARTING POINTS

* USING 1990 TO 2014

parameter income_shock_obs;
income_shock_obs(r)$pcgdp(r,"1990") = ((pcgdp(r,"2014") / pcgdp(r,"1990") -1 )/25) +1;
income_shock_obs("NIC") = income_shock_obs("HND");
income_shock_obs("ARE") = income_shock_obs("QAT");

* USING 2005 TO 2014
parameter income_shock_obs_2005;
income_shock_obs_2005(r)$pcgdp(r,"2005") = (pcgdp(r,"2014") / pcgdp(r,"2005") -1 ) +1;
income_shock_obs_2005("NIC") = income_shock_obs_2005("HND");
income_shock_obs_2005("ARE") = income_shock_obs_2005("QAT");



* PREDICTED GROWTH

* average 2018 to 2020 predicted growth rate
* computed in 'global economic prospects ..... .xlsx' excel file
* combination of World Bank and OECD forecasts

* https://data.oecd.org/gdp/real-gdp-forecast.htm#indicator-chart
* http://www.worldbank.org/en/publication/global-economic-prospects

parameter growthforecasts/
AUS	2.965
NZL	3.085
CHN	6.265
HKG	6.000
JPN	1.000
KOR	3.205
MNG	6.067
TWN	5.000
KHM	6.733
IDN	5.255
LAO	6.800
MYS	5.100
PHL	6.667
SGP	1.500
THA	3.900
VNM	6.633
BGD	6.900
IND	7.465
NPL	4.300
PAK	5.400
LKA	4.600
CAN	1.985
USA	2.414
MEX	2.531
ARG	2.670
BOL	3.633
BRA	2.793
CHL	3.396
COL	3.067
ECU	1.533
PRY	4.233
PER	3.700
URY	3.100
VEN	-8.433
CRI	3.554
GTM	3.233
HND	3.633
NIC	4.533
PAN	5.600
SLV	2.233
AUT	1.964
BEL	1.591
CYP	1.500
CZE	3.018
DNK	1.653
EST	3.246
FIN	2.181
FRA	1.771
DEU	1.791
GRC	2.762
HUN	3.440
IRL	3.040
ITA	1.029
LVA	3.490
LTU	2.730
LUX	3.335
MLT	1.500
NLD	2.666
POL	3.693
PRT	2.137
SVK	4.059
SVN	3.311
ESP	2.139
SWE	2.380
GBR	1.352
CHE	2.022
NOR	1.700
ALB	3.533
BGR	3.667
BLR	2.700
HRV	2.700
ROU	4.567
RUS	1.536
UKR	3.833
KAZ	3.267
KGZ	4.667
ARM	4.033
AZE	2.933
GEO	4.767
BHR	1.967
IRN	4.133
KWT	2.800
QAT	2.933
SAU	2.114
TUR	4.876
ARE	4.033
EGY	5.433
MAR	3.400
TUN	2.900
CMR	4.100
CIV	7.267
GHA	6.333
NGA	2.233
SEN	6.867
ETH	9.733
KEN	5.833
MDG	5.333
MWI	4.233
MUS	3.967
MOZ	3.433
TZA	6.800
UGA	6.000
ZMB	4.467
ZWE	3.500
BWA	3.367
NAM	2.267
ZAF	2.194
ISR	3.673
OMN	2.567/;

* NOTE: THESE AREN'T PER CAPITA
parameter income_shock_proj;
income_shock_proj(r) = growthforecasts(r) / 100 + 1;


* output the GDP shocks to compute correlations
chg_energy_co2("income shock - obs gdp 1990",r)  = income_shock_obs(r);
chg_energy_co2("income shock - obs gdp 2005",r)  = income_shock_obs_2005(r);
chg_energy_co2("income shock - proj gdp",r)  = income_shock_proj(r);

* for weighing, also output the initial emissions
chg_energy_co2("co2 - bmk emissions - direct ele",r) = content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r);
chg_energy_co2("co2 - bmk emissions - mrio",r) = content("co2","cons - mrio", "data", r) *content("co2","income", "data", r);



* then, compute approximations:
* first, compute approximated change in emissions if demand was HOMOTHETIC
* elasticity of world emissions to world income growth
*= delta emissions / emissions   / (delta income / income)
* = sum(emissions * shock-1) / sum(emissions)   / sum(income * shock-1) / sum(income)

chg_energy_co2("co2 - change - obs gdp 90 - H APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) *
  (income_shock_obs(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs(r)-1)) ;


chg_energy_co2("co2 - change - obs gdp 90 - H APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) *
  (income_shock_obs(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs(r)-1)) ;

*
* and now the Non-homothetic equivalent

chg_energy_co2("co2 - change - obs gdp 90 - NH APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons - direct ele",r) *
  (income_shock_obs(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs(r)-1)) ;


chg_energy_co2("co2 - change - obs gdp 90 - NH APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons mrio",r) *
  (income_shock_obs(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs(r)-1)) ;


* --------------------------------------------------------------

* SAME FOR 2005 AS START POINT:

chg_energy_co2("co2 - change - obs gdp 05 - H APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) *
  (income_shock_obs_2005(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs_2005(r)-1)) ;


chg_energy_co2("co2 - change - obs gdp 05 - H APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) *
  (income_shock_obs_2005(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs_2005(r)-1)) ;

*
* and now the Non-homothetic equivalent

chg_energy_co2("co2 - change - obs gdp 05 - NH APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons - direct ele",r) *
  (income_shock_obs_2005(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs_2005(r)-1)) ;


chg_energy_co2("co2 - change - obs gdp 05 - NH APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons mrio",r) *
  (income_shock_obs_2005(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_obs_2005(r)-1)) ;

* --------------------------------------------------------------

* SAME FOR PROJECTIONS



chg_energy_co2("co2 - change - obs gdp proj - H APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) *
  (income_shock_proj(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_proj(r)-1)) ;


chg_energy_co2("co2 - change - obs gdp proj - H APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) *
  (income_shock_proj(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_proj(r)-1)) ;

*
* and now the Non-homothetic equivalent

chg_energy_co2("co2 - change - obs gdp proj - NH APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons - direct ele",r) *
  (income_shock_proj(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_proj(r)-1)) ;


chg_energy_co2("co2 - change - obs gdp proj - NH APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons mrio",r) *
  (income_shock_proj(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_proj(r)-1)) ;

* --------------------------------------------------------------

* FINALLY, TO COMPARE, WITH UNIFORM GROWTh
parameter income_shock_uniform(r);

 income_shock_uniform(r) = 1.1;

* should be the same as computed elsewhere -- almost the case -- not quite

chg_energy_co2("co2 - change - uniform growth - H APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) *
  (income_shock_uniform(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_uniform(r)-1)) ;


chg_energy_co2("co2 - change - uniform growth - H APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) *
  (income_shock_uniform(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 

  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_uniform(r)-1)) ;

*
* and now the Non-homothetic equivalent

chg_energy_co2("co2 - change - uniform growth - NH APPROX direct ele","world") =

  sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons - direct ele",r) *
  (income_shock_uniform(r)-1) )/ 
 
 sum(r, content("co2","cons - direct ele", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_uniform(r)-1)) ;


chg_energy_co2("co2 - change - uniform growth - NH APPROX mrio","world") =

  sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r) * chg_energy_co2("co2 - elast - cons mrio",r) *
  (income_shock_uniform(r)-1) )/ 
 
 sum(r, content("co2","cons - mrio", "data", r) *content("co2","income", "data", r)) * 


  sum(r,content("co2","income", "data", r))   /
 
  sum(r,content("co2","income", "data", r) *
  (income_shock_uniform(r)-1)) ;

* --------------------------------------------------------------


* reduce size of output
option kill= trade;
option kill= trade_sh_counterfactual;
option kill= trade_sh;

* for revision, outputing with estimated theta
execute_unload 'results%SLASH%co2stats_%demandest%_%spec%_%prodspec%_fitteddem_%fitdem%_ESTTHETA.gdx'
execute_unload 'results%SLASH%simulationresults_tomerge%SLASH%%spec%_%prodspec%_fitteddem_%fitdem%_ESTTHETA.gdx', chg_energy_co2, content;
	
*execute_unload 'estimates/co2stats_%demandest%_%spec%_%prodspec%_fitteddem_%fitdem%_ESTTHETA_10pctshock.gdx'
*execute_unload 'estimates/chg_energy_co2_%demandest%/%spec%_%prodspec%_fitteddem_%fitdem%_ESTTHETA_10pctshock.gdx', chg_energy_co2;

*execute_unload 'estimates/co2stats_%demandest%_%spec%_%prodspec%_fitteddem_%fitdem%_ESTTHETA_obsgdp.gdx'
*execute_unload 'estimates/chg_energy_co2_%demandest%/%spec%_%prodspec%_fitteddem_%fitdem%_ESTTHETA_obsgdp.gdx', chg_energy_co2;

